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Abstract 

We describe a single step, second-order accurate Godunov scheme for ideal MHD 
based on combining the piecewise parabolic method (PPM) for performing spatial 
reconstruction, the corner transport upwind (CTU) method of Colella for multidi- 
mensional integration, and the constrained transport (CT) algorithm for preserving 
the divergence-free constraint on the magnetic field. We adopt the most compact 
form of CT, which requires the field be represented by area-averages at cell faces. 
We demonstrate that the fluxes of the area-averaged field used by CT can be made 
consistent with the fluxes of the volume-averaged field returned by a Riemann solver 
if they obey certain simple relationships. We use these relationships to derive new 
algorithms for constructing the CT fluxes at grid cell corners which reduce exactly 
to the equivalent one-dimensional solver for plane-parallel, grid-aligned flow. We 
show that the PPM reconstruction algorithm must include multidimensional terms 
for MHD, and we describe a number of important extensions that must be made 
to CTU in order for it to be used for MHD with CT. We present the results of a 
variety of test problems to demonstrate the method is accurate and robust. 
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1 Introduction 



In recent years a variety of numerical algorithms for multidimensional magne- 
tohydrodynamics (MHD) based on Godunov's method have been developed 
[3,9,10,16,23,25,28]. There are two important extensions to the basic hydro- 
dynamical algorithm that are required for MHD. The first is an extension 
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of the Riemann solver used to compute the fluxes of each conserved quan- 
tity to MHD; the second is a method by which the divergence free constraint 
V • B = is imposed upon the numerically evolved magnetic field. Of the two, 
the latter has emerged as the more difficult to achieve. 

The fact that it is important to ensure the numerically evolved field satisfies 
the divergence free constraint was first noted by Brackbill & Barnes [4]. They 
pointed out that the Lorentz force is not orthogonal to B if V • B 7^ 0, and 
that this could lead to incorrect dynamics in a well defined test problem. More 
recently Toth [32] has shown that in some circumstances it is possible to get 
the wrong jump conditions across MHD shocks if the constraint is not satisfied 
(this is also evident in the method of [16]). 

Currently, there are three methods by which the divergence free constraint 
is apphed in Godunov schemes. The first is to use a Hodge projection to 
clean the magnetic field of any divergence after each time step (e.g. [2,9,32]). 
The second is to extend the system of conservation laws with an evolutionary 
equation for the divergence designed to minimize the accumulation of error in 
any one location. Examples of this approach include the eight- wave scheme [26] 
and the GLM-MHD scheme [12]. Finally, the third is to design the difference 
equations for the magnetic field to explicitly conserve magnetic fiux, and so 
preserve the divergence free constraint. The latter method termed constrained 
transport (CT) by [15] has proved successful in other MHD algorithms [6,31] 
and is the method adopted here. 

The most compact CT difference formulae arc built upon area-averaged mag- 
netic field components located at the faces of a grid cell, rather than volume- 
averaged field components located at grid cell centers (CT algorithms built 
upon cell-centered fields have been developed in [32], however they require 
averaging over a stencil which is larger than that used to compute the fiuxes) . 
The need for a staggered grid is often thought of as a disadvantage of CT. In 
fact, however, it refiects one of the most attractive properties of CT: the funda- 
mental conserved property of the magnetic field in MHD is the magnetic fiux 
(which is an area- rather than volume-average) and by design, CT conserves 
the magnetic fiux (and therefore it preserves the divergence free character of 
the magnetic field) in an integral sense, over the smallest discretization scale, 
the grid cell size. 

We would like to combine CT with a finite volume shock capturing method. 
It may at first seem inconsistent to build a numerical algorithm which mixes a 
finite volume approach (which conserves integrals of volume-averaged values) 
with CT (which conserves integrals of area-averaged magnetic fiuxes at grid 
cell faces). In this paper we show that, provided the fluxes of the volume- and 
area-averaged fields obey certain simple relationships, the finite-volume and 
CT approaches can be made consistent (see also [23]). Most importantly, the 
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relationships wc derive allow us to develop algorithms for constructing fluxes 
of the face-centered field (located at grid cell corners in 2D) from the fluxes 
of the volume-averaged field computed by a Riemann solver (located at grid 
cell faces) . We demonstrate that multidimensional algorithms developed in this 
way reduce exactly to the equivalent one-dimensional solver for plane-parallel, 
grid-aligned flow. This is one difference between the methods developed here 
and previous implementations of CT in Godunov schemes [3,11,28]. 

We combine the CT algorithms developed in this paper with the piecewise 
parabolic method (PPM) [8] using Roe's linearization as the MHD Riemann 
solver [5] . An essential ingredient of PPM is a spatial reconstruction step to 
compute time-advanced estimates of the conserved variables at grid faces. 
In this paper we show that for MHD, this reconstruction step must include 
multidimensional terms in the induction equation (used to reconstruct the 
transverse components of the fleld). We argue that dimensionally split MHD 
algorithms cannot preserve the divergence free constraint between each one- 
dimensional update, therefore we adopt the unsplit corner transport upwind 
(CTU) algorithm of Colella [7] to develop a multidimensional algorithm. How- 
ever, there are a number of important extensions that must be made to CTU 
to make it suitable for MHD with CT. These include using a CT update for 
the magnetic field during the predict step, and inclusion of multidimensional 
terms along with the transverse flux gradients used to predict multidimen- 
sional fluxes. We describe these extensions in detail. 

The resulting two-dimensional MHD PPM algorithm uses a single step update, 
is second order accurate, and is fully conservative. Hence it is ideally suited for 
use on a statically or adaptively reflned mesh. Moreover, the two-dimensional 
algorithm reduces exactly to the base one-dimensional algorithm for planar, 
grid-aligned flows. Although this paper will only describe the combination of 
CT with CTU, there is no reason why the CT algorithm described here could 
not be combined with other unsplit methods. To simplify the discussion we 
conflne ourselves to a two-dimensional algorithm in this paper. 

The paper is organized as follows. In §2 we review the flnite volume and 
CT methods in order to demonstrate the relationships between the area- and 
volume-averaged magnetic flelds and their fluxes. In §3.3 we use these re- 
lationships to derive algorithms for constructing the fluxes of face-centered 
area-averaged flelds needed by CT from the fluxes of cell-centered volume- 
averaged flelds returned by a Riemann solver. We present a single step flrst- 
order Godunov method and use it to test different algorithms for constructing 
the fluxes used in CT. This method forms the flrst half of the second-order 
CTU + CT integration algorithm developed in section §4. In §5 we present 
a variety of tests which demonstrate the linear and nonlinear behavior of the 
scheme. Finally, in §6 we conclude. 
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2 Constrained Transport in Finite Volume Schemes 



The equations of ideal magnetohydrodynamics (MHD) can be written in con- 
servative form as 



_^ + V-(pv)=0 (1) 

^ + V- (pvv - BB) + VP* = (2) 
at 

— + V-(vB-Bv)=0 (3) 

— + V-((£; + P>-B(B.v))=0 (4) 



where p is the mass density, pv is the momentum density, B is the magnetic 
field, and E is the total energy density. The total pressure P* = P + (B ■ B)/2 
where P is the gas pressure, and the total energy density E is related to the 
internal energy density e via 

£; = e + p(v-v)/2 + (B-B)/2 . (5) 



Throughout this paper we will assume an ideal gas equation of state for which 
P = (7 — l)e, where 7 is the ratio of specific heats. Unless otherwise stated, 
we take 7 = 5/3. None of the main results described in this paper depend 
directly upon the equation of state. Note also that we have chosen a system 
of units in which the magnetic permeability /i — 1. 

In addition to the evolutionary conservation laws equations (1) through (4), 
the magnetic field must also obey the divergence free constraint, i.e. V • B = 0. 
It is of paramount importance that the numerically evolved field satisfy this 
constraint at all times, otherwise, for example, the system of equations for the 
conservative variables is inconsistent with the same system written in terms 
of the primitive variables, i.e. (p, v, B, P). 

The algorithm described in this paper is built upon finite volume (FV) meth- 
ods, in which the conserved variables are averaged over grid cell volumes. On 
the other hand, the CT method is built upon area- averaging of the magnetic 
field, leading to difference equations for the magnetic flux through the surfaces 
of grid cells. In the following subsections, we briefly review the FV and CT 
schemes in order to arrive at consistent relationships between the volume- and 
area-averaged magnetic field and their associated numerical fiuxes. 
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2.1 Finite Volume Method 



Consider a regular, two dimensional cartesian grid with grid cell {i,j) centered 
at (xi.yj) and of size {Sx,5y). The conservative system of equations for ideal 
MHD can be written in vector form as 







(6) 



where 
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is the vector of conserved variables and 



fx 



pVx 

Pv'x + P* -X 

pVxVy 

PVxVz 





Bl 



BxBy 

BxBz 



VxBz - B^Vz 



[{E + P*)vx-B,{B-v) ) 



) fy 



fWy 

fXVyVx - ByBx 

pvl + P*- Bl 

pVyVz 



VyBx 



ByBz 
ByVx 



VyBz - ByVz 
[{E + P*)Vy - By{B ■ V) ) 



,(8) 



are the flux vectors. Integrating over the volume of grid cell {i,j) and the time 
interval 5t — — and applying Gauss's theorem we obtain 
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the integral form of the evolution equation. The conserved quantities 

yi+5y/2 Xi+Sx/2 
yi-Sy/2 Xi-5x/2 

are averaged over the grid cell volume and the fluxes 

yi+Sy/2 

fS11/2,j I / ^^^^^ ^ ^)^^^^ (11) 

yi-Sy/2 
t"+i Xi+5x/2 

a;i-<5a;/2 

are averaged over the surface area of a grid cell face and the time interval St. 
Typically one approximates the flux integrals in equations (11) and (12) to 
some order of accuracy, while maintaining strict conservation by evolving the 
conserved quantities through equation (9). Note when written in this form, the 
components of the flux vectors arc non-zero for the transverse components of 
the magnetic field only, meaning that directionally split updates of the volume 
averaged field based on these fiuxes will not generally satisfy the divergence 
free constraint between directional sweeps. This suggests directionally split 
algorithms are inappropriate for MHD. 



2.2 Constrained Transport Method 



In the CT method, the integral form of the induction equation is based on 
area rather then volume averages. Starting from the differential form of the 
induction equation, 

— + VxS^O (13) 



where the electric field £ — — v x B in ideal MHD, one may integrate over the 
bounding surface of a grid cell and use Stoke's Theorem to obtain 

RH+l —fin .L^/'P^+VS _ cn+1/2 \ /.-v 

^x,i±l/2,j - ^x,i±l/2,j ^ r {^z,i±l/2,j-l/2 ^z,i±l/2,j+l/2j K^"^) 



5y 

^y,i,j±l/2 — ^y,i,j±l/2 \^z,i-l/2,j±l/2 '^z,i+l/2,j±l/2j K^^J 
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as the integral form of the evolution equation. The magnetic field components 

B:,±i/2,j ^j- I B,{x, ± e)dy (16) 

yi-5y/2 

Ki,3±i/2 = ^ / By{x, ± 5y/2, e)dx (17) 

Xi—5x/2 

are averaged over the grid cell bounding faces and 

0^/2j±i/2 = ll j ^z{x, ± 5x/2, y, ± 8y/2, t)dt (18) 

is averaged over the time interval 5t. Note the fundamental representation of 
the magnetic field is an area-average at cell faces. Although CT-like difference 
formulae are possible based on volume averaged fields at cell centers [32] , they 
preserve a discretization of the divergence on a different (larger) stencil than 
used to compute the fluxes. 

Just as in finite volume methods, one typically approximates the electric field 
(flux) integral in equation (18) to some order of accuracy and applies equations 
(14) and (15) to evolve the magnetic field components in time. Nevertheless, in 
a manner exactly analogous to the finite volume method, conservation of mag- 
netic flux is strictly enforced, implying that the net magnetic charge interior 
to a grid cell vanishes at time if it did so at time t". As such, the preser- 
vation of V • B = (in an integral sense) in the CT method is as fundamental 
as, e.g., the conservation of mass in a flnite volume method. Moreover, in the 
CT approach, the magnetic field components are averaged over the smallest 
dimensional volume necessary so as to transform the differential equation into 
its integral form. In this way one maintains the maximal "point-wise" infor- 
mation possible, thereby minimizing the dissipation inherent in the averaging 
process. 



2.3 Consistency of the CT & FV Methods 



To build a numerical scheme based on the CT method for the magnetic fiux, 
and a FV method for the remaining conserved quantities, it is very important 
that the surface and volume averaged magnetic field components (and their 
fluxes) be coupled in a consistent manner. One common approach [11,28,3] 
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(which wc also follow) is to define the volume-averaged magnetic field com- 
ponents at cell centers as equal to the average of the area-averaged values at 
cell faces, i.e. 



^y,i,j = 2 (^^'^'J-V2 + Ki,j+i/2) (20) 

which is sufficient for second order accuracy. However, as shown below this 
choice implies a specific relationship between the numerical fluxes for the in- 
duction equation as integrated in the CT and FV formulations. 

The expressions which describe the coupling of the fluxes in the CT & FV 
methods are a direct result of requiring that equations (19) and (20) also hold 
at time t""*"^. Subtracting equation (19) from an equivalent expression at time 
and substituting equations (9) and (14) for the time differences of the 
volume- and area-averaged magnetic field respectively we find 



2 fT^n+l/2 _ pn+1/2 \_}^{cn+l/2 _ cn+1/2 \ x 

^Bx • {^y,i,j-l/2 ^y,i,j+l/2) " ^ \^ z,i-l/2,j-\/2 ^ z,i-l/2,j+\/2) l^-^J 

1 /c.n+1/2 _ r.n+1/2 \ 

^ 2 V^,j+l/2J-l/2 ^z,i+l/2,jJrl/2) ■ \^^) 

where e^^ is a unit vector for the component of the flux vector. It follows 
that 

^ pn+l/2 _ 1 (cn+\/2 cn+l/2 \ 

^Bx • ^y,i,j±l/2 - 2 \f'z,i-\/2,j±l/2 ^ ^ z,i+l/2,j±l/2) l^'^J 



which is consistent with the observation that the flux averages used in the finite 
volume method (equation 12) are spatial averages of the electric fields used in 
the constrained transport method (equation 18). Repeating this analysis for 
By we find 

^ pn+l/2 _ ~1 |^c"+V2 ,cn+l/2 \ 

^By • ^x,i±l/2,j - 2 \^z,i±l/2,j-l/2 ^z,i±l/2,j+l/2) l^^J 

where is a unit vector for the By component of the fiux vector. Function- 
ally, equations (23) and (24) imply that one must replace the Godunov fluxes 
for the volume averaged x- and y-components of the magnetic field with the 
average of the corner centered regardless of the details of the CT algorithm 
used to compute the latter. Thus, equations (23) and (24) can be thought of 
as a corrector step that makes the predicted Godunov fiuxes given by the 
Riemann solver consistent with the CT fluxes. Clearly, the CT algorithm used 
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to compute the comer centered Ez will directly impact the accuracy and sta- 
bility of the underlying Godunov scheme; in §3.2 we address the problem of 
constructing CT algorithms and the properties they should possess. 



3 First Order CT Godunov Scheme 

In this section we will construct and test a single step, two dimensional, first or- 
der integration algorithm based upon the piecewise parabolic method (PPM). 
The simplicity of this algorithm allows us to develop and test two of the most 
important elements of this paper: the calculation of the interface states in 
MHD and the systematic construction of CT algorithms. Moreover, the re- 
sulting algorithm (apart from the update step) is essentially the first half of 
the CTU -|- CT integration algorithm described in §4. 

In section 3.1 we describe the calculation of the interface states in the PPM 
algorithm for the system of ideal MHD. This interface state calculation in- 
volves a characteristic evolution of a dimensionally split system. In this step 
we will find the appearance of truly multidimensional terms in the induction 
equation which are proportional to dB^/dx and dBy/dy. When these terms 
are included in the dimensionally split, linearized system of equations which 
are used to perform the characteristic evolution, they take the appearance 
of "source terms". We present two gedanken experiments which show that 
these multidimensional terms are essential to accurately predicting the time 
evolution in the interface states. 

In §3.2 we address the question of the consistency of the CT algorithm with 
the underlying finite volume integration algorithm. Repeating the arguments 
which lead to equations (23) and (24), we present an example of a CT algo- 
rithm which has insufficient dissipation and fails to reduce to the underlying 
integration algorithm for plane-parallel, grid-aligned fiows. We proceed to de- 
scribe a systematic approach to constructing a CT algorithm which can be 
applied with any approximate Riemann solver and reduces exactly to the un- 
derlying finite volume integration algorithm for plane-parallel, grid-aligned 
fiows. In §3.3 we present a numerical study comparing three, surprisingly sim- 
ple, CT algorithms each of which differs only in its dissipation for truly mul- 
tidimensional problems. 

3.1 Calculating the Interface States 

The algorithms presented in this paper are built upon the the piecewise 
parabofic method (PPM). For a thorough discussion of PPM or its finear 
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variant PLM and its implementation we refer the reader to the excellent de- 
scriptions in [8.29.24]. Roughly speaking the PPM algorithm can be broken 
down into three steps: spatial reconstruction, characteristic evolution, and flux 
evaluation. The purpose of these first two steps is to calculate a one-sided esti- 
mate of the time averaged state at the left- or right-hand sides of a particular 
grid cell interface. With these interface states in hand, the interface fiux may 
be calculated via either an exact or approximate Riemann solver. 

The calculation of the interface states in PPM is performed in primitive vari- 
ables, and is a one-dimensional algorithm. However, in the two dimensional 

{x, y) system of equations for ideal MHD there appear terms proportional to 
dB^/dx and dBy/dy which are not present in the truly one-dimensional sys- 
tem. In primitive variables, these terms only appear in the induction equation, 
which in component form is 

+ — {VyB, - ByV,) = (25) 

dB d 

-Qf + Q^ - B.vy) = (26) 

dB 

-gf + Q^ i^^B, - B.v,) + ^ {vyB, - Byv,) = . (27) 



Using the magnetic charge constraint (V • B = 0) these terms can be elimi- 
nated from equation (27), giving 



However, no such simplification can be made to equations (25) and (26). It is 
natural to ask just how important these terms are and what role they play in 
the evolution of the magnetic field. A few gedanken experiments quickly show 
that they are absolutely essential, and at times are the dominant term in the 
equation. 

Perhaps the most trivial example of a case in which these terms play an im- 
portant role is for stationary solutions. As a concrete example, consider a 
circularly polarized Alfven wave oriented at some oblique angle to the grid. 
One may always choose a reference frame in which the Alfven wave is station- 
ary. Because the wave is oriented oblique to the grid, dB^/dx and dBy/dy 
are non-zero throughout the domain (except at extrema). In addition, for a 
standing Alfven wave, the velocities are of the order of the Alfven speed. In 
this case, the term Vx{dBy/dy) in equation (25) is not a small term and in 
fact, it must exactly balance the remaining term d/dy{vyBx) — By{dvx/dy) 
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in order to hold the stationary solution. A similar situation also holds for 
equation (26). 



As a second example, consider the simple advection of a magnetic field loop 
in the {x, |/)-plane. Specifically, let (p, P, v) = a constant with v = v^i, = 
0, and a circular magnetic field loop in the {x, y)-plane of sufficiently weak 
strength that P — 2P/B'^ » 1. This problem is equivalent to the advection 
of a passive scalar, the 2;-component of the magnetic vector potential. In this 
case equations (25) and (26) are to a very good approximation given by 



dB^ dB., , , 

- v^—^ = (29) 



dt dy 



Hence it is clear that for this particular problem the term Vx{dBy/dy) is not 
only important, but completely controls the evolution of the x-component of 
the magnetic field. 



We conclude that if the calculation of the interface states for multidimensional 
ideal MHD includes a characteristic evolution step, it is necessary to include 
the infiuence of the inherently multidimensional terms in the induction equa- 
tion. Since PPM reconstruction includes a characteristic evolution step, we 
have found the following modifications necessary for multidimensional MHD. 



We will restrict the description to a single spatial grid cell index i and con- 
sider the reconstruction process in the x-direction. We begin by calculating 
the primitive state vector, = {p, v^, Vy, v^, B^, By, B^, P}i and such that 
Vi — {Vi, Bx^iji associated with g,, the vector of the cell averaged conserved 
variables. Next we apply the PPM algorithm to calculate the interface states 
of Vi where the characteristic evolution step is calculated by solving 
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where 
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(33) 



The matrix A is hnearized about the state and the source term a is taken 

to be a constant with the only non-zero term evaluated as Vy^i{Bx^i+i/2 — 
Bx,i-i/2) / ^x. Note that equation (31) includes all of the terms from equation 
(26). Denote the interface states calculated in this procedure as ^'^d 

^-1/2 where the superscripts (L, R) denote the left or right hand side of the 
interface to which they are adjacent. The final step is to define the primitive 
states V^^y^ = (i^i+i/2> Bx,i+i/2) and V^'^^/^ = (V'i-i/2> Bx,i-i/2)- Note a further 
significant advantage of using face-centered (staggered) fields: the interface 
states of the longitudinal component of the magnetic field do not need to be 
reconstructed, and therefore will be continuous. Moreover, since monotonicity 
constraints associated with reconstruction are not applied, extrema in the 
longitudinal component of B at interfaces will be preserved. 

The calculation of the y-interface states follows this same procedure, with the 
matrix A replaced by the equivalent one- dimensional wave matrix for the y- 
direction and the source term a containing a non-zero entry for Bx equal to 
Vx{dBy/dy). 
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3.2 Constrained Transport Algorithms 



In this section we take up the problem of constructing CT algorithms. In order 
to identify the properties of a suitable CT algorithm, it is particularly inter- 
esting to consider the limiting case of plane-parallel, grid-ahgned flows. In this 
hmit dBx/dx = —dBy/dy = so that there is no longer a difference between 
area and volume averaged magnetic fields. Moreover, if for example d/dx = 
then the correct solution to the CT algorithm is readily obtained via symmetry, 
e.g. £z,i+i/2,j+i/2 = {£z,i,j+i/2 + £z,i+\,j+i/2) /2. When a CT algorithm reduces 
to this, or an equivalent expression, for plane-parallel grid-aligned flows we de- 
scribe it as being consistent with the underlying integration algorithm, since 
in this case it will give the identical solution as the underlying integration al- 
gorithm applied to the equivalent one-dimensional problem. Furthermore, we 
seek to construct CT algorithms which are compatible with any approximate 
Riemann solver, e.g. [5,10,14,20]. Hence, they should only depend upon the 
electric field in the flux vector, not on the structure of the waves which result 
from the solution of the Riemann problem. 



3.2.1 Arithmetic Averaging 

Perhaps the simplest, and most often suggested, CT algorithm is based upon 
averaging the face centered electric fields obtained from the underlying inte- 
gration algorithm, i.e. choose ^^,i+i/2j+i/2 = ^^,i+i/2j+i/2 where 

^z,i+l/2J+l/2 = ^ {^■z,i+l/2,j + ^■z,i+\/2,j+l + ^z,jj+l/2 + ^■z, 1+1,3+1/2) ■ (34) 

Unfortunately, this CT algorithm is not consistent with the underlying inte- 
gration algorithm for plane-parallel, grid-aligned flows. This behavior is most 
easily understood when the underlying flnite volume integration algorithm is 
unsplit. 

Consider a plane-parallel, grid-aligned flow in which d/dx = 0. It follows that 
£z,i,j+i/2 = Sz,i+i,j+i/2 and Sz,i+i/2,j = £z,i,j- Inserting these expressions into 
equation (34) we find 

1 1 

£z,i+l/2,j+l/2 = ^ {^z,i,j + £z,i,j+l) + -^£z,i,j+l/2 ■ (35) 

Contrast this with the correct solution, which by the assumption of planar 
symmetry is simply Sz,i+i/2,j+i/2 = £z,i,j+i/2- In order to assess the impact 
this CT algorithm has on the integration algorithm as a whole, let the FV 
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numerical flux j j+1/2 be written as 

Fy,i,i+i/2 = 2 {fyi%j) + fyi%j+i) + A,j+i/2(gi,j - qi,j+i)) , (36) 

where is the viscosity-matrix [18,14]. Contracting this expression with 

a unit vector e^^ to extract the y-flux of (remembering that fy{Bx) — S^) 
we have the FV numerical electric field 

^z,i,j+i/2 = 2 + ^z,i,j+i) + 2^Ba:Dij+i/2{qi,j - Qij+i) ■ (37) 

Recall that, as discussed in §2.3, for FV + CT schemes this electric field is 
essentially a predictor value. To obtain the corrector value we begin by insert- 
ing the FV numerical electric field in equation (37) into the E CT algorithm 
in equation (35) giving 

Ez,i+i/2,j+i/2 = 2{^z,i,j + ^z,i,j+i) + ^e^^ Aj+i/2(?jj - Qi,j+i) ■ (38) 

Note also that by symmetry we have 4,i-i/2,i+i/2 = 4,j+i/2,i+i/2- Applying 
equation (23) we obtain the corrector value of the FV numerical electric field 
Ez,i,j+i/2 = ^2,i+i/2j+i/2 given by equation (38). Comparing equations (37) 
and (38) we find that the numerical viscosity is reduced by a factor of 2. 
Hence it's clear that with the arithmetic average CT algorithm S, the solution 
algorithm does not reduce to the underlying integration algorithm for plane- 
parallel, grid-aligned fiows and the stability of this approach is questionable. 
The failure of this simple procedure to reduce to the underlying integration 
algorithm for plane-parallel, grid-aligned fiows can be traced back to the lack 
of a directional bias in the averaging formula. 

The arithmetic average CT algorithm formed the basis of an algorithm pro- 
posed by Balsara and Spicer [3]. The need for the CT algorithm to have 
a directional biasing was well understood by these authors. In their paper, 
they presented two switches which serve as local, multidimensional sensors for 
magnetosonic shocks. The authors then applied weighting coefficients which 
impart a directional bias to the CT algorithm. 

However, the recognition that in equation (38) the viscous fiux contribution 
to the corner value for in the CT algorithm is simply too small by a factor 
of two, suggests that by doubling it we could recover the proper directional 
biasing. To that end we define 

^z,i+l/2,j+l/2 = 2£^,i+l/2j+l/2 
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(39) 



which can be written out exphcitly as 




(40) 



Repeating the arguments leading to equation (35) using the CT algorithm 
defined by equation (40) one finds that the resulting scheme reduces to the 
underlying finite volume integration algorithm for plane-parallel, grid-aligned 
flows. It is not clear, however, from the ad-hoc construction described here 
how well such an algorithm will behave for truly multidimensional flows. 

In the following section we describe a systematic approach to constructing a 
CT algorithm which by design reverts to the underlying flnite volume inte- 
gration algorithm for plane-parallel, grid-aligned flows. Two entirely new CT 
algorithms will be constructed based upon different approximations. We will 
also find that the CT algorithm described by equation (40) can be understood 
as a limiting case of one of the CT algorithms constructed in the next section. 
Despite the simphcity of the CT algorithms constructed in the following sec- 
tion, we will see in §3.3 that they behave surprisingly well on multidimensional 
tests. 

3.2.2 Systematic Construction of CT Algorithms 

The approach described here is based upon the observation that a CT algo- 
rithm can be thought of as the inverse of the consistency relations given by 
equations (23) and (24). In this sense, the CT algorithm can be thought of 
as a reconstruction, or integration procedure. Note that the interface fluxes 
which are calculated as part of the base integration algorithm are midpoint 
values, centered spatially on the grid cell face, and averaged temporally over 
the time step. This suggests that we consider the CT algorithm, which cal- 
culates a time averaged value of at the grid cell corner, to be a spatial 
integration procedure. For example, given a face centered value £z,i+i/2,j we 
seek an estimate of (9£^/5y)i+i/2,j+i/4 giving one value for 



Clearly in two dimensions one may integrate from any one of the four nearest 
face centers to the corner and generally the resulting values for Sz,i+i/2,j+i/2 




(41) 
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will differ. In the CT algorithms presented here, we will use the arithmetic 
average of these four values giving 



^z,i+l/2,j+l/2 ^ 



+ 



- (^£z,i+l/2,j + ^z,i+l/2,j+l + £z,i,j+l/2 + ^z,i+l,j+l/2^ 

Sy I (dsA _ (ds^\ 

J i+l/2J+l/4 ) i+l/2J+3/4y 



° V \ "-^ / i+l/4j+l/2 \ / i+3/4j+l/2y 

The construction of this CT algorithm is completed by specifying a way to 
calculate the derivatives of £z on the grid cell face. 

To calculate {d£z/dx) and (dSz/dy) at grid cell faces, we propose to use an 
approximate solution for the evolution equations for (dB^/dx) and (dBy/dy). 
At a |/-interface we differentiate the induction equation for giving 

d_ fdBA ^ d_ (d£A _ ^ 
dt \ dx j dy \dx j 



Similarly, at an x-interface we differentiate the induction equation for By giv- 
ing 



Since these expression are still in conservation form it suggests that we may 
calculate an interface value for {d£z/dx) at y-interfaces and {d£z/dy) at x- 
intcrfaces using for example an HLL or Lax-Friedrichs flux. To evaluate these 
fluxes, we need estimates for the gradients of (dSz/dx), {dSz/dy), [dB^/dx) 
and (dBy/dy) on either side of the interface. 

For the single step, CT Godunov algorithm which we are considering in this 
section, we calculate these derivatives as follows. For {dB^/ dx) we difference 
the interface and cell center values giving 




For (dSz/dx) we difference the face centered Sz,i+i/2,j which comes directly 
from the Riemann solver and the cell center value £z,i,j evaluated in the cell 
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center state g^^- giving 



TZ {^z,i+i/2,j - E;,,i,j) ■ (46) 



dx I . . 6x 



The values for {dBy/dy) and {dSz/dy) are given by analogous expressions. 
Pursuing the Lax-Priedrichs estimate with a maximum wave speed a we find 



J i+l/4,i+l/2 



r- {Ez,i+l/2,j - Ez,i,j + £z,i+l/2,j+l - ^z,i,j+l) 



+ {Bx,i+l/2,j - Bx,i,j - -Ba:,i+l/2j+l + -Ba;,ij+l) (47) 



and 



( ) - TT. {^z,i,j+l/2 - Ez,i,j + Ez,i+l,j+l/2 - ^z,i+l,j) 

J i+i/2,i+i/4 

+ ^ {By,i+i,j+i/2 - By,i+ij - + . (48) 

Repeating this procedure for the two remaining gradients and inserting the 
results into equation (42) we obtain 



^z,i+l/2,j+l/2 — 2 {^z,i,j+l/2 + Ez,i+l,j+l/2 + £z,i+l/2,j + £z,i+l/2,j+l) 

^ ^z,i+l,j ~l~ Ez,i,j+1 ~l~ ^2:,i+lJ+l) 

+ {Bx,i+l/2,j - Bx,i,j - Bx,i+l/2,j+l + -Ba;,iJ+l) 

+ (-6a;,i+l/2j ~ Bx^i+lJ — Bx^i+l/2,j+l + B^^i+lj+lj 

+ (-Sy,i+lj+l/2 - By^i+ij - By^ij+i/2 + -By.ij) 

+ g- (-Bj,,i+ij+i/2 - - By^ij+i/2 + By^iJ+l^ . (49) 

One may readily show that for plane-parallel, grid-aligned flows, this CT al- 
gorithm will properly recover the associated one-dimensional solution for the 
underlying integration algorithm. Hereafter, we refer to equation (49) as the 
£^ CT algorithm. 

It is particularly interesting to note that the a = limit of equation (49) 
gives equation (40). Hence we may now understand equation (40) as being 
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equivalent to the integration and averaging procedure described here with the 
approximation 

I dy J ,+1/2 j+i/4 ^\9yJ ,,,+1/4 ^\dy) ,+i,,.+i/4 ' 



Clearly this is not an upwinded approximation, suggesting that we should find 
some level of oscillations present in using this CT algorithm for multidimen- 
sional flows. However, since the dissipation arising from the terms proportional 
to a in the CT algorithm are only important for truly multidimensional 
fiows, it is not clear if their neglect will have a substantive impact on the first 
order integration algorithm which we are considering here. For that reason, 
henceforth wc will refer to the a = limit of the £^ CT algorithm as S°, and 
include it in the tests in the following section. 

We now have two CT algorithms: the £° algorithm given by equation (40), 
and the £^ algorithm given by equation (49). As our final CT algorithm, we 

note that for the special case of advection, {dSz/dy) at an x-interface should 
be selected in an upwind fashion according to the contact mode. As such, 
we suggest that upwinding {d£z/dy) at x-interfaces (and similarly {d£z/dx) 
at y-interfaces) according to the contact mode may be sufficient to lead to a 
stable, non-oscillatory integration algorithm. Specifically, we choose 

{dSJdy)ij+i/4 for ^;x,j+i/2j > 

{d£Jdy)i+ij+i/4 for v^^i+i/2,j < (51) 

+ {dSJ dy)i+ij+i/i^ otherwise . 



, dy 



i+l/2j+l/4 



Note that this simply depends upon the sign of the mass fiux, not the details 
of the solution of the Riemann problem at the interface and therefore can be 
applied with any approximate or exact Riemann solver. An analogous expres- 
sion holds for the remaining three interface gradients of Sz- We will refer to 
the CT algorithm which results from combining this approximation for the 
gradients of £z with equation (42) as the £^ CT algorithm. By design this CT 
algorithm reduces to the underlying integration algorithm for plane-parallel 
grid-aligned flows and is properly upwinded in a multidimensional sense for 
the simple case of magnetic field advection. 

In this section we've presented a simple approach to constructing a CT al- 
gorithm which reduces exactly to the base integration algorithm for plane- 
parallel, grid-aligned flows. By design, the CT algorithms constructed here 
differ only in their numerical viscosity for multidimensional problems. It should 
be noted that the approach presented here can readily be incorporated into 
other numerical schemes, such as wave propagation algorithms [19]. The CT 
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algorithms described here can also be applied to integration algorithms based 
upon spatial reconstruction to the grid cell corners as well [1]. While this might 
seem surprising following this presentation, note that the factors {5x, 6y) can- 
cel in the £°, Ef) CT algorithms. For such integration algorithms, the cell 
center value £z,i,j should be replaced with the value of £z calculated in the 
reconstructed fluid state at the grid cell corner, e.g. qi+i/2,j+i/2- With this 
choice, the CT algorithms presented here will reduce to the base integration 
algorithms for plane-parallel, grid-aligned flows. 

We note that another CT algorithm with the properties that it reduces to the 
base integration algorithm for plane-parallel, grid-aligned flows has recently 
been presented and tested elsewhere [23] . In particular, the authors of that pa- 
per present a general framework for combining CT and Godunov-type schemes 
and two specific implementations for their positive and central-type schemes. 
A direct comparison between their approach and ours is somewhat complex 
in the general case. In the specific case of a first order Godunov scheme, one 
can show (using equations 41 - 47 in [23]) that their CT algorithm is identical 
to the £1 CT algorithm constructed here, although this is not immediately 
obvious from the description of their framework. For the more complex CT 
algorithms developed in our paper, it is likely that they too can be cast in 
the framework described by Londrillo & Del Zanna [23] , although we have not 
attempted to do so. 



3.3 First Order CT Godunov Tests 



The integration algorithm utilized in this section is easily assembled from the 
elements described in the preceding sections. Starting at time we calculate 
the X- and y-interface states as described in section 3.1 using the PPM algo- 
rithm. Next we use an approximate Riemann solver to calculate a flux at each 
grid cell interface. In the tests presented here we use a Roe linearization [5]. 
Finally, we apply one of the three {£^,£°,£^) CT algorithms. The resulting 
integration algorithm is flrst order accurate and subject to a restrictive CFL 
stability limit. In the tests presented in this section we use a time step 



where A™""^ indicates the fastest wave mode speed in the x- or y-direction. 




(52) 
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3.3.1 Field Loop Advection 



The first problem we consider is the advection of a weak magnetic field loop. 
The computational domain extends from —l<x<l, and —0.5 < y < 0.5, 
is resolved on a 2N x N grid and has periodic boundary conditions on both 
X- and y-boundarics. For the tests presented in this section we take = 64. 
The mass density p = 1 and the gas pressure P = 1. The velocity components 
Vx = Vocos{9), Vy = fosin(6'), Vz = where cos(6') = 2/\/5 and sin(6') = l/\/5. 
In the diffusion tests we set Vq — while in the advection tests we set Vq — y/5 
so that by t = 1 the field loop will have been advected around the grid one 
complete orbit along the grid diagonal. The 2;-component of the magnetic 
field Bz = while the in plane components and By are initialized from the 
2;-component of the magnetic vector potential where 



A., = I 



Ao(R -r)iorr<R 

^ ' - (53) 

for r > i? 



where Aq = 10"^ R = 0.3 and r = ^Jx^ + y\ Thus for r < i?, /? = IPjB'^ = 
2 X 10^ and the magnetic field is essentially a passive scalar. 

In the first test we consider the diffusion of the field loop. In figure 1 we 
compare grey-scale images of the magnetic pressure (5^ + -^^) at t = to the 
evolved results at t = 2 using the three Ef) CT algorithms. Clearly the 

E'^ CT algorithm leads to an unacceptable amount of diffusion compared to 
the other two, as evidenced by the emergence of a hole at the center caused by 
reconnection. The (£^°,£^) CT algorithms lead to essentially identical results 
and a very small amount of diffusion. Apparently the additional dissipation 
included in the CT algorithm is not necessary for stability in this test. 

Next we consider the advection of the magnetic field loop. Unlike in the sta- 
tionary field loop test, here the evolved results appear quite similar for the E'^ 
and E'z CT algorithms. In figure 2 we present grey-scale images of the mag- 
netic pressure at t = 0.19 for the three CT algorithms. At this 
time, the E'^ and E^ CT algorithms give quite similar results. In contrast, the 
E° CT algorithm appears to have insufficient dissipation leading to an oscilla- 
tory solution. These observations are consistent with the comments in §3.2.2 
regarding the upwinding of the gradients of Ez at the interfaces. In figure 3 
we present grey-scale images of the magnetic pressure at t = 2 for the three 
CT algorithms. By this time the oscillations present using the E^ 
CT algorithm have come to dominate the solution. The results from the £° 
and El CT algorithms continue to remain quite similar, implying that the 
dissipation of the first order scheme is comparable to the dissipation in the 
CT algorithm. Note also the similarity to the magnetic pressure image in 
figure 1 for the E^ CT algorithm. 
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Fig. 1. Grey-scale images of the magnetic pressure (i?^ + ^^) at t = (top left) and 
at i = 2 for a stationary medium (vq = 0) using the (top right), the £° (bottom 
left) and the £^ (bottom right) CT algorithm. 




Fig. 2. Grey-scale images of the magnetic pressure (B^ + By) at t = 0.19 for an 
advected field loop {vq = \/5) using the £" (top left), £° (top right) and £^ (bottom) 
CT algorithm. 

Prom the results of these tests we are led to conclude that of the three al- 
gorithms, the CT algorithm is preferable. The CT algorithm leads to 
stable, yet diffusive results for stationary problems. The S° CT algorithm 
appears to have insufficient dissipation for advection problems leading to os- 
cillatory results. One might naturally wonder, however, if these results are 
biased to favor the £^ CT algorithm by design. In the following section we 
present additional tests of these CT algorithms where wave modes other than 
the contact mode play an important role in the solution. We note in passing 
that the source terms described in §3.1 are absolutely essential to obtain the 
results presented here. If they had been omitted, the field loop disintegrates 
in oscillations before completing a fraction of an orbital period. 
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Fig. 3. Grey-scale images of the magnetic pressure (B'^ + By) at t = 2 for an advected 
field loop {vq = -\/5) using the £^ (top left), £° (top right) and £^ (bottom) CT 
algorithm. 

3.3.2 Circularly Polarized Alfven Wave 

In a recent paper Toth [32] described a test problem involving the evolution 
of traveling and standing circularly polarized Alfven waves in a periodic do- 
main. This test problem is interesting from the point of view that the initial 
conditions are nonhnear solutions to the equations of ideal MHD. Unfortu- 
nately, their efficacy as a discriminating test for multidimensional MHD codes 
has been hindered slightly [25,23] by the fact that they are susceptible to a 
parametric instability [17,13]. Nevertheless, we have found this to be a useful 
test and find no indication of instability for the parameters adopted here. 

The initial conditions we utilize here are slightly different than in the original 
description [32]. The computational domain extends from < a; < y^, and 
< ?/ < \/5/2, is resolved on a 2N x N grid and has periodic boundary 
conditions on both x- and |/-boundarics. For the tests presented in this section 
we take N — 8. The Alfven wave propagates at an angle 9 — tan~^(2) 63.4° 
with respect to the x-axis and has a wavelength A = 1. The mass density p — 1 
and the gas pressure P = 0.1. The velocity and magnetic field components are 
most easily described in a rotated coordinate system 



xi = ,T cos 6' + y sin 6' (54) 
a;2 = — a; sin 6* -|- 7/ cos ^ (55) 
X3 = z (56) 

such that the Alfven wave propagates along the Xi axis. The magnetic field 
components Bi = 1, B2 = 0.1 sin(27ra;i), and B3 = 0.1 cos(27rxi). The velocity 
components Vi = (0, 1) for traveling or standing Alfven waves respectively, 
V2 — 0.1 sin(27rxi), and vs — 0.1 cos(27ra;i). With this set of initial conditions 
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and vi = the Alfven wave will travel a distance of one wavelength A in a 
time t — 1. 

To better illustrate the geometry of this problem, a high resolution image of 
the out of plane component of the magnetic field, is presented in figure 
4. Also included in this figure is an image of for the resolution tested 
(A^ = 8) in order to emphasize that the coarsest resolution for this wave is in 
the y-direction (with only 8 grid points per wavelength), and that there are 
essentially two complete wavelengths across the grid diagonal. We have found 
that low resolution tests such as presented below are much more informative, 
since differences between algorithms are generally largest in this case. Note 
that, just as in the field loop problem, the in-plane components of the magnetic 
field {Bx, By) are initialized via the 2;-component of the appropriate magnetic 
vector potential. 




Fig. 4. Plot of Bz at the initial time for the circularly polarized Alfven wave problem 
at high resolution, = 128, (left) and at the resolution plotted in figure 5 (right). 

In Toth's analysis [32] he found that the errors in the solution were dominated 
by errors in the transverse magnetic field B2 and velocity f 2 in our notation. 
He presented line plots of B2 versus a; as a function of resolution and numerical 
scheme. In figure 5 we present analogous line plots of B2 versus Xi for the case 
of a standing and traveling waves including the initial conditions at t — and 
the solutions at t — 5 for the three CT algorithms under study. It is worth 
noting that in these plots we have calculated B2 using the cell center magnetic 
fields and have included the data for every grid cell in the calculation. Owing 
to the angle of the wave with respect to the grid, there are many data points 
which fall on a single xi position. The lack of scatter in the data shows that 
the Alfven wave retains its planar structure extremely well for the duration of 
the calculation. 

Comparing the plots of B2 versus Xi for the three CT algorithms, we generally 
find that the £° or £^ CT algorithm give nearly identical results, while the £^ 
CT algorithm is more dissipative. In the case of the standing wave solution, 
the dissipation rate in the CT algorithm is approximately twice that of 
the other two. In the traveling wave case, the difference in the dissipation 
rate is much less indicating that the dissipation in the first order integration 
algorithm is comparable. 
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standing Alfven Wave 



Traveling Alfven Wave 




Fig. 5. Plot of B2 versus Xi for the standing (left) and traveling (right) circularly 
polarized Alfven wave problem at t = and t = 5 for the three {£^,S°,S^) CT 
algorithms. 

From these tests, and many additional tests not included here, we conclude 
that the £^ CT algorithm has the best dissipation properties. The £^ CT 
algorithm gives stable, yet diffusive results while the £° CT algorithm appears 
to have insufficient dissipation, leading to oscillatory results for advection 
problems. These observations suggest that the CT algorithm with a set 
equal to some estimate of the magnitude of the local gas velocity could also 
lead to a non-oscillatory CT algorithm. This idea is however untested, nor is 
it clear that it would result in an algorithm which is superior to the CT 
algorithm. For the remainder of this paper we will use the £^ CT algorithm 
in test problems. 



4 Second Order CTU + CT Godunov Scheme 



The corner transport upwind (CTU) method was developed by Colella as an 
unsplit, two dimensional sequel to the piecewise parabolic method (PPM) of 
Colella and Woodward for Euler's equations. The algorithm is second order 
accurate and degenerates to the base PPM algorithm for plane-parallel, grid- 
aligned flows. Formally, the algorithm can be described in just a few steps. 

First, one calculates left and right interface states at each grid cell face using 
the 1-dimensional algorithm from the base PPM scheme. Using the notation 

adopted in §3.1, let these be denoted by (gf+i/2j, 9fl*i/2,j^ ^iS+i/^)- 
For MHD, the PPM interpolation used to construct these states must include 
all the multidimensional terms identified in §3.1. At each interface one solves 
the Riemann problem associated with these interface states and computes the 
fluxes {F*^^_^_^/2,j^ Py,i,j+i/2)- Next, one updates the interface states to the 1/2 
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time step 



QM-l/2,j = Q4+l/2,j ^ 2dy ('^?'*»+l'i+V2 ~ ^y,i+l,j-l/2) (58) 
?i\i+l/2 = 1i,j+l/2 ~ 2^ (-^^Vl/2j ~ -^x*i-l/2j) (59) 
(li^j+1/2 = ?i\j+l/2 ~ 2 5x ("^^'^+1/2 J+1 ~ -^x,i-l/2J+l) (60) 



Solving the Riemann problem associated with the four updated interface states 
i%+i/2ji ^ij+1/2) '^^'^ obtains second order accurate fluxes which can be used 
to update qij via the standard finite volume integration relation, equation (9). 

Unfortunately, the CTU method as just described is incomplete for MHD 
when using constrained transport. One reason is that equations (57) - (60) 
fail to preserve the V ■ B = condition. The marriage of CTU with CT for 
MHD requires a modification of equations (57) - (60) for updating the interface 
states, and an additional CT integration step for updating q'^j from time 
tot"+^ 



4-1 Updating the Interface States 



There are two modifications to the CTU method required for MHD using CT. 
The first modification, required by constrained transport, is that the fiuxes 
i'^xi+i/2j' Pyij+1/2) niust be integrated from face center, to the grid cell cor- 
ner. Hence from the interface centered flux i^^j_|_i/2 j obtain two corner cen- 
tered fluxes {,Fx^+i/2,j+i/2i Fx^+i/2,j-i/2) where the superscript L, R indicates 
that the fiuxes have been integrated from face center to the grid cell corner 
from either the left (— y) or right (+?/) side of the grid cell corner. The labehng 
of the y-fiuxes (Fj;*+i/2j+i/2 ' -^3^1-1/2^+1/2) obtained from i^*^jj_,_i/2 follows an 
analogous convention. The second modification, originating from differences 
in the form of the equations of ideal MHD when written in primitive, or con- 
servative variables, is the addition of source terms. The modified form of the 
update relations in equations (57) - (60) for advancing the interface states to 
the 1/2 time step can be formally written as 

1i+i/2,j = <li+i/2,j ~ 2 5y (■^^»*+V2j+i/2 ~ -^^^^+1/2^-1/2) + -^^x,i,j (61) 
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Qi+l/2,j - 1i+l/2,j 2 Sy V 2/'»+V2j+l/2 ^y,i+l/2,j-l/2) + g ^''''+^'^ 
?i\i+l/2 = ~ 2^ (-^^i+l/2j+l/2 - Px,i-l/2,j+l/2) + -^Sy^i,j (63) 

(li!j+l/2 = ~2'5x i^Si+l/2,j+l/2 - F^,i-l/2,j+l/2) + -^Sy,ij+i . (64) 



The procedure for integrating the fluxes from face center to grid cell corner, 
and the need for the source terms will now be described in turn. 



4.1.1 Integrating the Fluxes to the Grid Cell Corner 

Recall that fx {By) = —£z and fy{Bx) = £z and let e^^ and denote the 
unit vectors for the and By components of the flux vector. Furthermore, 
let ^:,i,j+i/2 = • F; . .^^^^ and £li+if2,j = -es, • ^;,i+i/2j- This allows us to 
define 

^y,i+l/2,j+l/2 = Py,i,j+l/2 + {^z,i+l/2,j+l/2 ~ ^z,i,j+l/2)(^B^ (65) 
^^i+l/2,j+l/2 = Fy,i+l,j+l/2 + (^2*j+l/2j+l/2 " ^2,1+1,^+1/2)^5^ (66) 
-^x,i+l/2,j+l/2 = -^a;,i+l/2,j ~ i^z,i+l/2,j+l/2 ~ ^z,i+l/2,j)(^By (67) 
-^i^«*+l/2,j+l/2 = -^x,i+l/2,j+l ~ (^2,i+l/2,j+l/2 ~ ^z,i+l/2,j+l)^By ■ (68) 

From a practical point of view, integrating the fluxes from face center to grid 
cell corner in this fashion can be thought of as simply stating that the normal 
components of the magnetic field at grid cell interfaces is advanced to time 
fn+i/2 ^-g^ ^ jjj^ggpa^i Namely, 

-^r,t±l/2,j = -S",i±l/2,j + 2 5y (^*,»±V2j-l/2 - ^r,i±l/2,j+l/2) (69) 
^yXj±l/2 = ^^,i,j±l/2 ~ 2 5x (^2*»-V2,i±l/2 - ^2,i+l/2,j±l/2) • (^O) 



The last part of this integration procedure which requires description is the 
CT algorithm used to calculate the corner centered emf ^2^+1/2^+1/2- To ac- 
complish this, note that the fluxes (-^^,+1/2^) ^yij+1/2) equivalent to the 
fluxes used in the single step integration algorithm tested in §3.3. As such, the 
calculation of the corner centered emf 4+1/2^+1/2 be accomplished with 
any of the {8^,8°, S^) CT algorithms described in §3.2. For the tests prob- 
lems presented in the following sections we will use the 8^ CT algorithm. We 
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note that following this procedure, the magnetic fields satisfy the V • B = 
condition at time i""*"^/^. 



4.1.2 Interface State MHD Source Terms 

The source terms present in equations (61) - (64) follow from the recognition 
that if they are set to zero, the updated interface states are not formally 
advanced by 5t/2 for MHD. The basic reason for this discrepancy lies in the 
fact that the interface states as described in §3.1 (and typically implemented 
in PPM) are calculated in primitive variables. As a result. 



2 dx 



(71) 



and similarly for the other interface states. To correct this situation, we define 

/ n \ 



Q . . — 







-"x,i+l/2,j -"a:,i-l/2j 

8x 



(72) 



and 



x,i,3 
^z,i,3 



z,i,3 



\ z,h3 z,i,3 I 



^y,i,j+l/2 ^y,i,j-l/2 



(73) 
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With this choice, the interface states as updated by equations (61) - (64) 
include all of the necessary terms so as to be advanced to time t""*"^/^. 



Note that the choice to include the term Vz{dBx/dx) for (and the associ- 
ated energy source term) at the x-interfaces in this step, rather than includ- 
ing it when calculating the interface states as described in §3.1 has a very 
important consequence; it prevents an erroneous field growth of in certain 
circumstances. To elucidate this situation, consider a magnetic field loop in 
the {x, |/)-plane advected with a uniform v — and set B^ — initially. 
With the exception of extrema, dB^/dx and dBy/dy are non-zero through- 
out the field loop. However, since Vz is uniform the magnetic field Bz should 
remain equal to zero. If the term VzdB,j./dx is included when calculating the 
x-interface states in §3.1 they would contain a non-zero Bz- Owing to the co- 
herent structure of the in-plane field loop, the the values of Bz in the interface 
states will also have a coherent structure. Upon updating the interface states 
with the transverse fiux gradients, the growth of Bz is diminished, however it 
is not canceled identically. The net result is an unphysical growth of a coher- 
ent Bz which eventually infiuences the in-plane dynamics. The choice of source 
terms described in section §3.1 and this section maintains Bz to the level of 
roundoff error with an incoherent structure. Hence, the algorithm presented 
here accurately captures the balance of the terms proportional to dB^jdx 
and dBy/dy^ as described in §3.1, in both the predictor and corrector steps 
for calculating the interface state values of Bz- 

4-2 The Constrained Transport Update Algorithm 

After having updated the interface states to time t^+^Z^ via equations (61) - 
(64), the interface fiux calculation is repeated giving rise to the second order 
accurate fiuxes (-^"^1/2 j' ^yt]+i/2)- ^^'^ CTU algorithm, this set of fiuxes 
is used to evolve q^- to time However, in order to evolve the magnetic 

fields via constrained transport, we must extend the CT algorithms described 
in §3.2. Requiring that the algorithm reduce to the base integration algorithm 
for plane-parallel, grid-aligned flows we flnd that we simply need to advance 
the electric field gradient calculation to the half time step, i.e. equation (46) 
is replaced with 




(74) 



The electric field ^"^jj is the cell center value advanced by 5t/2, that is 



.n+1/2 



= V. 



n+1 2 nn+1/2 



— V. 



n+1/2 j^n+1/2 



(75) 



28 



where to be consistent with the integration scheme, the cell center magnetic 
fields are given by 



pn+l/2 _ 1 /pn+l/2 nn+1/2 \ 

^y,i,j - 2 \^y,i,j-l/2 + ^y,i,j+l/2) \") 

where the field components on the right hand side of these equations are 
equal to the normal components of the magnetic field in the interface states 

(li±i/2j^ 1i 'j±i/2)- "^^^ density, x- and y- momenta needed to compute the 
velocity components in equation (75) are advanced by dt/2 using 

QiJ ^ = +2^ {K,i-l/2,j - Fx,i+l/2,j) + 2^ (-^y*ij-l/2 ~ Py,i,j+l/2) (78) 

where the fiuxes are those calculated in the first step of the CTU algorithm. 
4-3 Summary 

The following steps summarize the CTU + CT algorithm for MHD: 

Calculate the x- and y-interface states {qf*i/2 j, ?i+i/2 i' ^i^7*+i/2> lfj+1/2) 



(2 

(3 
(4 
(5 

(6 

(7 
(8 



j+l/2j' yj+l/2,j> yi,j+l/2> Hi,j+l/2) 

using the PPM algorithm, and the multidimensional source terms as de- 
scribed by equations (31) - (33) in §3.1. 

Calculate the x- and y-interface fiuxes (-^^^i+i/2j, -^^,1^+1/2) associated 

with the interface states (?^i/2 j' ^1*3+1/2) ^ Riemann solver. 
Using the CT algorithm described in equations (42) and (51) integrate 
the face centered fluxes to the grid cell corner as described in §4.1.1. 
Compute the the four updated interface states (5j^f/2 j' li j+1/2) equa- 
tions (61) - (64) with the source terms detailed in §4.1.2. 
Compute the x- and y-interface fiuxes {F^t+i%j^ ^ytj+1/2) associated 
with the interface states (?j^^/2i' ^ij+1/2) ^ Riemann solver. 

Compute the grid cell corner centered electric field i+i/2 i+1/2 ^^ing the 
CT algorithm described in equations (42) and (51) advanced to time 

^n+1/2 described in §4.2. 

Advance the surface averaged normal components of the magnetic field 
from time to t'^'^^ using equations (14) and (15). 
Advance the remaining volume averaged conserved quantities from time 
to t'"'^^ using equation (9). 



This completes the description of the algorithm. It is second order accurate, 
unsplit, and preserves the V • B = constraint throughout the time step. In 
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the following section we apply this CTU + CT scheme to a variety of test 
problems. 



5 Tests 

In this section we present results obtained with the CTU + CT integration 
algorithm just described. Throughout these tests we use the CT algorithm. 

5.1 Field Loop Advection 

The advection of a magnetic field loop discussed in §3.3.1 was instructive for 
assessing the dissipation in the different CT algorithms. In this subsection we 
present the results obtained for this problem with the second-order CTU + 
CT algorithm. The grid resolution, and initial conditions are equivalent to 
those used in §3.3.1. 

In figure 6 we present grey-scale images of at times t = and t = 2. 
Comparing these figures we find that the majority of the field dissipation has 
occurred at the center and boundaries of the field loop, where the current 
density is initially singular. A more quantitative measure of the magnetic field 
dissipation rate is given by the time evolution of the volume average of 
as shown in figure 7. We find that the measured values (denoted by symbols) 
is well described by a power law (solid line) of the form B"^ = A{1 — {t/r)"-) 
with A = 3.463 x 10-^ r = 10.614 x 10^ and a = 0.2914. 

Another important indicator of the properties of the integration algorithm 
is the geometry of the magnetic field lines. Note that since the CT method 
evolves the interface magnetic flux (preserving V • B = 0) one may readily 
integrate to find the ^-component of the magnetic vector potential. The mag- 
netic field lines presented in figure 8 are obtained by contouring A^. The same 
values of A^ are used for the contours in both the t = and the t = 2 images. 
By t = 2 the inner most field line has dissipated. It is quite pleasing, however, 
to note that the CTU -|- CT algorithm preserves the circular shape of the 
magnetic field lines, even at this low resolution. 

5.2 Circularly Polarized Alfven Wave 

The test problem involving the propagation of circularly polarized Alfven 
waves at an oblique angle to the grid was described in §3.3.2. In this sub- 
section we present a resolution study for both standing and traveling Alfven 
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Fig. 6. Grey-scale images of the magnetic pressure (B^ + By) at i = (left) and 
t = 2 (right) using the CTU + CT integration algorithm. 
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Fig. 7. Plot of the volume averaged magnetic energy density B^ as a function of 
time. The solid line is a power law curve fit to the data points denoted by symbols. 




Fig. 8. Magnetic field lines at t = (left) and t = 2 (right) using the CTU + CT 
integration algorithm. 

waves. The initial conditions are equivalent to those used in §3.3.2 only with 
N = {4, 8, 16, 32}. 

As a diagnostic of the solution accuracy, we plot the in-plane component of 
the magnetic field, B2, perpendicular to the wave propagation direction, xi, 
in figure 9. These plots are constructed using the cell center components of 
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the magnetic field, and each grid cell is inchided in the plots. Hence, the lack 
of scatter demonstrates that the solutions retain their planar symmetry quite 
well. Figure 9 includes the solutions at time t = 5 with N = {4, 8, 16, 32} 
for both standing and traveling waves. For comparison, we also include the 
initial conditions for the = 64 case. We find that these results compare well 
against other published calculations [32,25] and we find no indication of the 
parametric instability for the parameters adopted here. 



Standing Alfven Wave Traveling Alfven Wave 




Fig. 9. Plot of B2 versus Xi at t = 5 for the standing (left) and traveling (right) 
circularly polarized Alfven wave problem. For comparison, the initial conditions at 
t = for the iV = 64 case is also included. 



5.3 Rotated Shock Tube Problem 

The solution to the one dimensional Riemann problem has long been used 
as a test of numerical algorithms [30,21]. Solving the same problem in a two 
dimensional domain with the initially planar discontinuity rotated by some 
angle with respect the the grid can also be a robust test of the integration 
algorithm. In addition to the usual questions, one is also interested in how 
well the planar symmetry is preserved for fiows which are oblique to the grid. 

For MHD this implies a particularly stringent condition on the component of 
the magnetic field in the direction of the initial discontinuity normal. Using 
the coordinate transformations in equations (54) - (56), let the initial discon- 
tinuity lie in the plane xl = constant. Then for MHD, the solution to the 
one-dimensional Riemann problem should have Bi = constant, which requires 
a balance between the x- and y-gradients of such that d£z/dx2 = 0. Toth 
[32] has recently shown that in some cases, schemes which do not preserve 
the V • B = condition can result in a solution in which Bi contains a jump 
across a shock; see also [16]. 

In the trivial case, rotated shock tube problems are initialized with the shock 
tube discontinuity oriented at a 45 degree angle with respect to the grid, i.e. 
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with a coordinate rotation angle 6 = ta.n~^(5x/Sy). Examples of test cal- 
culations performed using this configuration can be found in [11,28,16]. We 
have run a variety of shock tube problems with this configuration (with both 
5x — 5y and 5x ^ 6y) and find that in all of our tests, the parallel compo- 
nent of the magnetic field, Bi, remains equal to a constant with variations 
which are of the order of roundoff error. We assert that this is a result of the 
symmetry of the initial conditions with respect to the grid. 

A non-trivial configuration with a coordinate rotation angle of ^ = tan~^(2) 
63.4° and Sx — Sy was recently suggested by Toth [32] and has been adopted 
elsewhere [23,9] as well. This problem is more challenging because at the dis- 
crete, grid-scale level the initial conditions contain variations along the plane 
of the initial shock tube discontinuity. Moreover, the symmetry in this config- 
uration is such that Qij = gj+2j_i which is outside of the integration stencil 
for most integration algorithms. This is especially true in the neighborhood 
of shocks where most integration algorithms drop to first order. Typically, it 
is in the neighborhood of shocks where one find oscillations, or in some cases 
even jumps, in Bi. 

We choose to simulate rotated shock tube problems on a grid of x Ny grid 
cells with Sx = Sy and the shock tube discontinuity oriented along the grid 
diagonal. Let C equal the greatest common divisor of (A^^:, Ny) and define 
= Nx/C and Vy = Ny/C. With this configuration, the coordinate rotation 
angle 9 = ta.n^^{rx/ry) and the symmetry is such that Qij = qi^r^,j-ry- Note 
that this computational grid can also be described as containing C x C "macro- 
cells" each of which is Tx x Ty grid cells in size. We have run a variety of shock 
tube problems with (r^, r-y) = (2, 1), (3, 2), (5, 4), etc. and in all cases find 
results which are mutually consistent. 

In the interest of presenting solutions which can be compared to previously 
published results [23,9,32] we will now focus on the {r^, Vy) = (2, 1) case with 
Nx = 256 and Ny = 128. The particular problem studied has a left state given 
by = (1, 10, 0, 0, 5/V471, 5/v/47r, 0, 20) and a right state given by = 
(1, -10, 0, 0, 5/v^,0, 1) where V = (p, vi, V2, v^, Bi, B2, B^, P). 

Among other places, the one dimensional solution to this Riemann problem 
can be found in figure la of [27]. 

In figure 10 we present line plots of the parallel component of the magnetic 
field Bi versus the parallel coordinate Xi. These line plots include every point 
in the computational domain, hence the lack of scatter indicates that the 
solution retains the planar structure quite well. The first line plot, labeled 
"grid cell" , is constructed using the cell center magnetic field components. We 
find oscillations in Bi which are roughly 10% of Bi, with the largest oscillations 
occurring at the fast-mode shocks and weaker oscillations at the left and right 
propagating slow-mode rarefaction and shock respectively. We note that the 
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CT algorithm does not reduce these oscillations further when compared to 
the CT algorithm. Hence the oscillations in Bi are not a result of insufficient 
dissipation in the CT algorithm. The second fine plot, labeled "macro-cell", 
is constructed by first conservatively averaging the solution onto a grid of 
128 X 128 "macro-cells" before computing the macro-cell center component of 
Bi. The variations in Bi when averaged onto a macro-cell arc of the order of 
roundoff error. Note that we obtain the same result for other rotation angles, 
e.g. with {r^, Vy) = (3, 2), (5, 4), etc. We conclude that the oscillations in Bi 
versus Xi are a simple consequence of the fact that on the scale of grid cells, 
the discretized solution contains variations in the X2-direction. Upon averaging 
the solution onto the macro-cells, this variation is eliminated, and we recover 
the condition Bi = a, constant. Note that this also suggests that if one wishes 
to eliminate the oscillations in Bi it would require a viscosity with a stencil 
whose size is at least as large as the macro-cell. 

The recovery of i?i = a constant upon averaging the solution onto a grid of 
macro-cells is clearly consistent with magnetic flux conservation and plane 
parallel symmetry, yet it does not appear to be a trivial result. For example, 
it is clear that schemes which generate a jump in Bi [32,16] can not recover 
this result. 



grid cell 

— macro-cell 

1.5 - 



1.2 I ' ^ ' ^ ' ^ ' ^ ' 1 

0.2 0.4 0.6 0.8 1 

Fig. 10. Plot of B\ versus x\ for the (r^;, Vy) = (2, 1) case at time t = 0.08 using 
the grid ceU centered B and the macro-cell centered B. The data for the macro-cell 
centered B has been offset vertically by 0.01 to for clarity. 



5.4 Linear Wave Convergence 



In this subsection we show that the CTU -|- CT integration algorithm con- 
verges with second order accuracy for linear amplitude waves. The computa- 
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tional domain extends from < x < 2/^/5, and < y < 1/^/5, is resolved 
on a 2N x grid and has periodic boundary conditions on both x- and y- 
boundaries. The hnear wave propagates at an angle 6 = tan^-'^(2) ~ 63.4° with 
respect to the a;-axis and has a wavelength A = 2/5. Using the coordinate ro- 
tation described by equations (54) - (56), the initial conserved variable state 
vector is given by 

= q + eRk cos(27rxi) (79) 



where q is the mean background state, s = 10~^ is the wave amplitude, and Rk 
is the right eigenvector in conserved variables for wave mode k (calculated in 
the state q). In order to enable others to perform the same tests presented here 
and compare the results in a quantitative manner, we include the numerical 
values for the right eigenvectors in the appendix. As in previous 2D calcula- 
tions, the in-plane components of the magnetic field {B^, By) are initialized 
via the 2;-component of the magnetic vector potential. 

The mean background state q is selected so that the wave speeds are well 
separated and there are no inherent symmetries in the magnetic field orienta- 
tion. It is most convenient to describe it in terms of the associated primitive 
variables and in the rotated coordinate system given by equations (54) - (56). 
The density p = I and gas pressure P = I/7 = 3/5. The velocity component 
parallel to the wave propagation direction, vi — 1 for the entropy mode test 
and ■Ui = for all other wave modes. The transverse velocity components 
V2 = = 0. The magnetic field components Bi = 1, B2 = \/2, and B^ = 1/2. 
With this choice, the slow mode speed Cg = 1/2, the Alfven speed Ca — 1, and 
the fast mode speed c/ = 2 in the wave propagation direction. 

The error in the solution is calculated after propagating the wave for a dis- 
tance equal to 1 wavelength. Hence, the initial state is evolved for a time 
t = \/c where c is the speed of the wave mode under consideration. For each 
component k of the conserved variable vector q we calculate the LI error with 
respect to the initial conditions 

= ^ EE (80) 



by summing over all grid cells We use the cell center components of the 

in-plane magnetic field components (S^, By) in computing this error. In figure 
11 we plot the norm of this error vector 



\M = jy:w (81) 
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for the fast, Alfven, slow and entropy modes. This plot shows that the solution 
for each wave mode converges with at least second order accuracy. The order 
of convergence for each wave mode, obtained by a power law fit to the errors, 
is indicated in the legend of figure 11. We note in passing that if the interface 
state reconstruction algorithm is performed using piecewise linear interpola- 
tion, instead of piecewise quadratic, the error is proportional to N""^ for all 
wave modes and the amplitude is increased slightly. 

Linear Wave Convergence 
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Entropy Mode -O (2.11) 
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100 



Fig. 11. Linear wave convergence of fast, Alfven, slow and entropy modes using 
the CTU + CT integration algorithm. The symbols denote the calculated LI error 
norm. The lines are power law curve fits to the data and the order of convergence 
for each wave mode is indicated in the legend. 



5.5 Current Sheet 



In this subsection we present a problem which is particularly sensitive to the 
numerical dissipation and demonstrates of the robustness of the integration 
algorithm. The computational domain extends from < a; < 2, and < |/ < 2, 
is resolved on an 256 x 256 grid and has periodic boundary conditions on both 
X- and y-boundaries. The density p — 1 and the magnetic field components 



— B, 



and 



B,, 



Bo if0<x<l/2 
-Bo if l/2<x< 3/2 
Bo if3/2<x<2 



(82) 
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where Bq — 1. Hence there arc initially two current sheets in the computa- 
tional domain and the characteristic Alfven speed = Sq/^/P = 1- The gas 
pressure P = 0.1 such that /3 = 2P/ Bq = 0.2 and the dynamics are initially 
magnetically dominated. The ratio of the Alfven speed to the sound speed 
-Bo/VtP ~ 2.45, hence magnetically driven dynamics are supersonic. The 
initial velocity components Vx — i'osin(7ry) with vq = Vy = Vz = For 
Vo/ca the ensuing dynamics are well characterized by linear Alfven waves. 
For the values selected here this is approximately true at early times, until 
magnetic reconnection and nonlinear effects come to influence the dynamics. 

One aspect of this problem which is of particular interest is the magnetic re- 
connection since it is a direct measure of the numerical resistivity. In figure 12 
we present the time evolution of the magnetic field lines. From the magnetic 
field geometry at time t — 0.5 we see that, as one should expect, the numerical 
resistivity is a function of the magnetic field orientation with respect to the 
grid: magnetic fields dissipate, and reconnect preferentially where the mag- 
netic field orientation is oblique to the grid. Hence, we find the largest change 
in the magnetic field structure at the nodal points of the transverse velocity. 
As reconnection takes place, the magnetic energy is converted into thermal 
energy (on time scales of the integration time step) which in turn drives both 
compressional, and Alfvcnic waves. These waves interact seeding more recon- 
nection events. By time t = 1 a series of magnetic islands have developed 
along the current sheets. These islands are free to move parallel to the local 
magnetic field direction and by time t — l.b some of the magnetic islands have 
propagated toward the velocity anti-nodes and merged. This process of island 
formation, translation, and merging continues until there are two magnetic 
field islands along each current sheet located approximately at the velocity 
anti-nodes. 

This problem is also interesting in that it uses a very simple set of initial 
conditions to test the "robustness" of the integration algorithm. The nonlinear 
dynamics which result from this problem lead to strong compressions and 
rarefactions. It is important to maintain the divergence free constraint as the 
topology of the field changes during reconnection. By either increasing vq or 
decreasing P (and therefore (5) the dynamics become increasingly difficult for 
the integration algorithm to solve. We have found it a very useful test to 
discriminate between algorithms. 



5.6 MHD Blast Wave 

As our final test problem we consider the explosion of a centrally over pressur- 
ized region into a low pressure, low (3 ambient medium. This problem has been 
studied by a number of authors [33,3,22] and we've chosen to use the param- 
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Fig. 12. Time evolution of the magnetic field lines using the CTU + CT integration 
algorithm. Time increases from left to right and top to bottom in normal reading 
order. The contour levels of Az which are plotted is uniform over the sequence of 
images at times t = (0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0). 

eters given by [22]. The computational domain extends from —0.5 < x < 0.5 
and —0.5 < y < 0.5. The density p = 1, the velocity v = 0, and the magnetic 
field components Bx = By = 10/\/2 and B^ = 0. Within a circle of radius 
R = 0.125 about the origin the gas pressure P = 100 and /3 = 2P/B'^ = 2. 
Outside of this circle, the gas pressure P — 1 and P — 2 x 10"^. 

The solution to this problem at time t = 0.2, using a 200 x 200 grid, is 
presented in figure 13. The density image shows two dense shells of gas which 
propagate parallel to the magnetic field. The outer surface of these shells is 
a slow-mode shock and the inner surface is the contact surface separating 
the gas initially inside and outside of the boundary surface. The maximum 
compression of the gas is 3.3 indicating that the slow mode shock is quite 
strong, in agreement with what one might expect from the ratio of the gas 
pressures, Pin/Pamb = 100. In the direction orthogonal to the magnetic field, 
the magnetic pressure is the dominant player in the dynamics, yet from the 
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field lines we see that there is only a moderate change in the geometry of the 
field. 



The solution presented in figure 13 demonstrates that the algorithm presented 
in this paper is both stable and accurate for low-/3 plasma problems involving 
strong MHD shock waves. The solutions also preserve the initial symmetry of 
the flow exceptionally well despite the orientation of the magnetic field. We 
find no indication of grid related artifacts in the solution. 




Fig. 13. Linearly scaled grey-scale images and the magnetic field lines of the evolved 
state (time=0.2) for the MHD blast wave problem. The density (top left) ranges 
from 0.192 (white) - 3.31 (black). The gas pressure {top right) ranges from 1.0 
(white) - 32.1 (black). The magnetic energy density {bottom left) ranges from 23.5 
(white) - 77.7 (black). 



39 



6 Conclusion 



In this paper we have demonstrated that the method of constrained transport 
can be combined with finite volume integration algorithms in a self consistent 
manner. Consistency, however, implies that the electric fields used in the CT 
update step and those used for evolving the volume average magnetic fields 
are coupled. This coupling has direct consequences for the stability and ac- 
curacy of the integration algorithm. We have presented a general approach 
to constructing CT algorithms, constructing and testing three. Each of these 
CT algorithms contained the novel property that for planar, grid-aligned fiows 
the solution would recover the ID solution obtained with the underlying in- 
tegration algorithm. These CT algorithms differed only in their dissipation 
properties for truly multidimensional fiows. Through numerical experiments 
we have shown that the CT algorithm is well behaved leading to stable, 
non-oscillatory solutions. We have also noted how this algorithm can readily 
be combined with other unsplit integration algorithms such as central schemes, 
or wave propagation methods. 

We have shown that if the PPM integration algorithm is used for ideal MHD, 
terms proportional to dBx/dx and dBy/dy, which in primitive variables are 
present only in the induction equation, must be included in the calculation of 
the "interface states". If these terms are neglected in the calculation of the 
interface states, the integration algorithm is oscillatory for the simple case 
of field loop advection. We also presented two simple gedanken experiments 
to demonstrate why this result should be expected. A simple approach for 
including these source terms in the calculation of the interface states is adopted 
and it has been demonstrated to be accurate and stable for use with both the 
single step and two step (CTU -|- CT) integration algorithm. 

Another result of this paper is the extension of the CTU integration algorithm 
for ideal MHD based upon the CT integration algorithm. We showed that since 
the interface states are calculated using primitive variables, the standard CTU 
procedure for updating the interface states to the 1/2 time step is missing 
terms which are proportional to dBx/ dx at x-interfaces and similarly for the 
y-interface states. These terms must be included as an additional set of "source 
terms" so that the interface states are formally advanced to the 1/2 time step. 
We also described how the CT algorithms developed in this paper can be 
combined with the CTU integration algorithm so as to maintain V ■ B = 
throughout the integration time step. 

The CTU + CT integration algorithm presented in this paper for ideal MHD 
has been thoroughly tested and some representative solutions have been in- 
cluded here. This algorithm combines the strong stability and shock capturing 
characteristics of Godunov methods with the magnetic fiux conservation ob- 
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tained via the CT method. The integration algorithm is conservative, uses a 
single step update algorithm, and is second order accurate on smooth solu- 
tions. These characteristics make it ideally suited for use on either a statically 
or adaptively refined mesh. We note also that "physical" source terms such as 
an external gravitational field, or Coriolis terms accounting for a rotating ref- 
erence frame can be readily incorporated into this integration algorithm. The 
resulting algorithm also retains the desirable properties noted above, such as 
recovering the ID solution for plane parallel grid-aligned flows. These details 
will be described in a later paper. 

Lastly, while the algorithm presented in this paper has focused solely on the 
two-dimensional case, the results presented here can be extended to three 
dimensions. This extension principally involves modifications to the "source 
terms" involved in the calculation and update of the interface states with 
transverse flux gradients. The actual details of this extension are beyond the 
scope of this paper and will be presented elsewhere. 
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A Linecir Wave Right Eigenvectors 



In order to enable others to perform the linear wave convergence test presented 
in section 5.4 and compare their results in a quantitative manner, we include 
the numerical values for the right eigenvectors here. In the rotated coordinate 
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system described by equations (54) - (56) the conserved variable vector 
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The right eigenvectors (labeled according to their propagation velocity) are 
given by 
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